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Abstract 

The cumulative effects of weak resonant and secular perturbations by the 
major planets produce chaotic behavior of asteroids on long timescales. Dy- 
namical chaos is the dominant loss mechanism for asteroids with diameters 
> 10 km in the current asteroid belt. In a numerical analysis of the long 
term evolution of test particles in the main asteroid belt region, we find that 
the dynamical loss history of test particles from this region is well described 
with a logarithmic decay law. In our simulations the loss rate function that 
is established at t ~ 1 My persists with little deviation to at least t = 4 Gy. 
Our study indicates that the asteroid belt region has experienced a signif- 
icant amount of depletion due to this dynamical erosion — having lost as 
much as ~ 50% of the large asteroids — since 1 My after the establishment 
of the current dynamical structure of the asteroid belt. Because the dynam- 
ical depletion of asteroids from the main belt is approximately logarithmic, 
an equal amount of depletion occurred in the time interval 10-200 My as 
in 0.2-4 Gy, roughly ~ 30% of the current number of large asteroids in the 
main belt over each interval. We find that asteroids escaping from the main 
belt due to dynamical chaos have an Earth impact probability of ~ 0.3%. 
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Our model suggests that the rate of impacts from large asteroids has declined 
by a factor of 3 over the last 3 Gy, and that the present-day impact flux of 
D > 10 km objects on the terrestrial planets is roughly an order of magni- 
tude less than estimates currently in use in crater chronologies and impact 

hazard risk assessments. 
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1. Introduction 



The main asteroid belt spans the ~ 2-4 AU heliocentric distance zone that 
is sparsely populated with rocky planetesimal debris. Strong mean motion 
resonances with Jupiter in several locations in the main bel t cause asteroid s 



to follow chaotic orbits and be removed from the main belt (jWisdom 



19871) 



These regions are therefore emptied of asteroi ds over the 



tem, forming the well-known Kirkwood gaps ( IKirkwoodl . 118671 ). In addition 



age o 



the solar sys- 



to the well known low-order mean motion resonances with Jupiter that form 
the Kirkwood gaps, there are numerous weak resonances t hat cause long term 



orbital chaos and transport asteroids out of the main belt (IMorbidelli and Nesvorny 



19991 ). A very powerful secular resonance that occurs where the pericenter 
precession rate of an asteroid is nearly the same as that of one of the solar 



system's eigen frequencies, the z/r secu 



the main belt (jWilliams and Faulkner! . Il98ll ). 



ar reso nance, lies at the inner edge of 



The many resonances found throughout the main asteroid belt are largely 
responsible for maintaining the Near Earth Asteroid (NEA) population. Non- 
gravitational forces, such as the Yarkovsky effect, cause asteroids to drift 
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in semimajor axis into chaotic resonances whence they can be lost from 



the main belt (lOpikl. 



1998 



Bottke et al. 



1951 



Vokrouhlicky and Farinellal . 



2000 



Farinella et al. 



20001 ). The Yarkovsky effect is size- dependant, and there- 
fore smaller asteroids are more mobile and are lost from the main belt more 
readily than larger ones. Asteroids with D < 10 km have also undergone ap 



preciable col 



1996 



Cheng 



isiona 



2004 



evolution over the a ge of the solar system (IGeissler et al. 



Bottke et al. 



20051 ), and collisiona. 



ject fragments into chaotic resonances (IGladman et al. 



events can also in- 



19971 ). These pro- 



cesses (coUisional fragmentation and semimajor axis drift followed by injec- 
tion into resonances) have contributed to a quasi steady-state flux of small 
asteroids (-D < 10 km) into the terrestrial planet region and are responsi- 
ble for del ivering the rnajority of ter re strial planet impac tors over the last 



3.5 Gy flBottke et al. 



200C 



2002a 



Strom et al. 



20051). 



In contrast, most members of the population of D > 30 km asteroids 
have existed relatively unchanged, both physically and in orbital properties, 
since the time when the current dynamical architecture of the main asteroid 
belt was established: the Yarkovsky drift is negligble and the mean coUi- 
sional breakup time is > 4 Gy for D > 30 km asteroids; asteroids with 
diameters between ~ 10-30 km have been moderately altered by coUisional 
and non-gravitational effects. However, as we show in the present work, the 
population of large asteroids {D > 10-30 km) is also subject to weak chaotic 
evolution and escape from the main belt on gigayear timescales. By means 
of numerical simulations, we computed the loss history of large asteroids in 
the main belt. We also computed the cumulative impacts of large asteroids 
on the terrestrial planets over the last ~ 3 Gy. 
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Our present study is additionally motivated by the need to understand 
better the origin of the present dynamical structure of the main asteroid belt. 
The orbital distribution of large asteroids that exist today was determi i ied by 



dynamical processes in the early solar system. iMinton and Malhotral (120091 ) 
showed that large asteroids do not uniformly fill regions of the main belt 
that are stable over the age of the solar system. The Jupiter-facing bound- 
aries of some Kirkwood gaps are more depleted than the Sunward bound- 
aries, and the inner asteroid belt is also more depleted than a model asteroid 
belt in which only gravitational perturbations arising from the planets in 
their cu rrent orbits have sculpted a n initially uniform distribution of as- 



teroids. 



Minton and Malhotral (120091 ) showed that the pattern of depletion 



observed in the main asteroid belt is consistent with the effects of resonance 



early in solar systern history (iFernandez and Ip 



Hahn and Malhotra 



1999 



Gomes et al 



. 1984: 


Malhotra . 


1993. 


1995 



20051 ). and that this event was the 



last major dynamical depletion event experienced by the main belt. The last 
major dynamical depletion event in the main asteroid belt likely coincided 
with the so-called Late Heavy Bombardment (LHB) ~ 3.9 G y ago as indi- 



cated by the crater record of the inner planets and the Moon (jStrom et al. 



20051). 



Knowledge about the distribution of the asteroids in the main belt just 
after that depletion event may help constrain models of that event. Quanti- 
fying the dynamical loss rates from the asteroid belt also help us understand 
the history of large impacts on the terrestrial planets. Motivated by these 
considerations, in this paper we explore the dynamical erosion of the main 
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asteroid belt, which is the dominant mechanism by which large asteroids have 
been lost during the post-LHB history of the solar system. 

We have performed n-body simulations of large numbers of test particles 
in the main belt region for long periods of time (4 Gy and 1.1 Gy). We 
have derived an empirical functional form for the population decay and for 
the dynamical loss rate of large main belt asteroids. Finally, we discuss the 
implications of our results for the history of large asteroidal impacts on the 
terrestrial planets. 



2. Numerical simulations 

Our long term orbit integrations of the solar system used a parallelized im- 
plementation of a second-orde r mixed variable symp l ectic mapping known as 



the W isdom- Holman Method (jWisdom and Holman 



1991 



Saha and Tremaine 



19921 ). where only the massless test particles are parallelized and the mas- 
sive planets are integrated in every computing node. Our model included 
the Sun and the planets Mars, Jupiter, Saturn, Uranus, and Neptune. All 
masses and initial conditions were taken from the JPL Horizons servicqj on 
July 21, 2008. The masses of Mercury, Venus, Earth, and the Moon were 
added to the mass of the Sun. 

Test particle asteroids were given eccentricity and inclination distribu- 
tions similar to the observed main belt, but a uniform distribution in semi- 
major axis. The initial eccentricity distribution of the test particles was 
modeled as a Gaussian with the peak at = 0.15 and a standard deviation 



^see http://ssd.jpl.nasa.gov/7horizons 
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of 0.07, a lower cutoff at zero, and an upper cutoff above any value which 
would lead to either a Mars or Jupiter-crossing orbit, whichever was smaller. 
The initial inclination distribution was modeled as a Gaussian with the peak 
/i = 8.5° and standard deviation of 7°, and a lower cutoff at 0°. The other 
initial orbital elements (longitude of ascending node, longitude of perihelion, 
and mean anomaly) were uniformly distributed. The eccentricity and inclina- 
tion distributions of the adopted initial conditions and those of the observed 
asteroids of absolute magnitude H < 10.8 are shown in Fig. [H 

Mars was the only terrestrial planet integrated in our simulations. De- 
spite its small mass. Mars has a significant effect on the dynamics of the inner 
asteroid belt due to numer ous weak resonances, including t hree-body Jupiter- 



m 



19991 ). Two simulations 



Mars-asteroid resonances flMorbidelli and Nesvornyl . 
were performed: Sim 1, with 5760 test particles integrated for 4 Gy, and 
Sim 2, with 115200 test particles integrated for 1.1 Gy. In each of these sim- 
ulations an integration step size of 0.1 y was used. Particles were considered 
lost if they approached within a Hill radius of a planet, or if they crossed 
either an inner boundary at 1 AU or an outer boundary at 100 AU. 

We define time t = as the epoch when the current dynamical architec- 
ture of the main asteroid belt and the major planets was established. What 
we mean by this is the time at whi ch any primordial m ass depletion and exci- 



tation has already taken place (see 



O'Brien et al. 



migration of giant plane t s has finished (IFernandez and id. 



1993 



Strom et al. 



2005 



Gomes et al. 



20071). and 



any early orbital 



2005 



1984 



Ma: 



Minton and Malhotra 



hotral . 



20091). 



At this epoch the main belt would have already had its eccentricity and in- 
clination distributions excited by some primordial process, and its semimajor 
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Figure 1: Eccentricity and inclination distributions. a),b) The distributions of the 931 
observed bright {H < 10.8) asteroids in the main belt that are not members of coUisional 
famil ies, from the AstDys online data service ( Knezevic and Milani . ! 



20061 ). c),d) The initial e and i distributions for Sim 1 (the 5760 particle simulation). e),f) 



The initial e and i distributions for Sim 2 (t^e 115200 particle simulation). 



axis distribution shaped by early planet migration. Therefore, the e and i 
distributions at t = likely resembled those of the present-day asteroid belt, 
although subsequent long-term evolution likely altered them somewhat from 
their primordial state. Current understanding of planetary system formation 
suggests that the epoch prior to when we define t = could have been several 
million to several hundred million years subsequent to the formation of the 
first solids in the protoplanet ary disk; the first solids have radiometrically 



determined ages of 4.567 Gy (IRussell et al. 



20061). 



3. Main asteroid belt population evolution 

The loss history of particles from Sim 1 and Sim 2 are shown in Fig. [2j 
The loss histories are nearly indistinguishable over the 1.1 Gy length of Sim 2. 
The loss history appears to go through two phases. The first phase, lasting 
until ~ 1 My, is characterized by a rapid loss of particles from highly unstable 
regions, such as the major Kirkwood gaps and the uq secular resonance. The 
slope of the loss rate on a log-linear scale changes rapidly between 0.3-1 My 
until the second phase is reached, which lasts from 1 My until at least the end 
of Sim 1 at 4 Gy. The slope of the loss rate on a log-linear scale continues to 
change during the second phase, but only over much longer timescales and 
by a much smaller amount than during the first phase. 

The particle removal times and particle fates (whether they become inward- 
going Mars-crossers, or outward-going Jupiter-crossers) for Sim 2 are both 
shown in Fig. [31 We find that the particles that are lost during the initial 
1 My (the red to light-green points) are generally those with high initial 
eccentricity, particles from the resonance (appearing as a curving yel- 
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low band in the semimajor axis vs. inclination plot), and particles in the 



strongly chaotic mean motion resonances wi th Jupiter (the Kir 



Michtchenko et al. 



swood gaps). 
f 2009L how- 



These maps are similar to those produced by 
ever their maps were coded by spectral number (more chaotic orbits having 
a larger spectral number) using 4.2 My integrations. The apparent rapid 
change in slope around 10^-10^ years is likely due mostly to the emptying of 
asteroids from the i/g resonance region (see also the upper right-hand panel 
of Fig. [3D. Also, in the outer asteroid belt the sheer proximity to Jupiter 
and the resulting strong short-term perturbations cause particles to be lost 
very rapidly. Many of these regions may never have accumulated asteroids, 
and therefore the loss of these particles represents a numerical artifact in our 
simulation due to over-filling the model asteroid belt with test particles. For 
instance, if the asteroid belt formed with the giant planets in their current 
positions, those regions would always have been unstable to asteroids, so 
none could have formed there. However, regions of the asteroid belt that are 
currently highly unstable may not always have been. Models of early solar 
system history indicate large changes to the orbital properties of the giant 



plane ts (IFernandez and Ip . 



1984 



Hahn and Malhotra 



1999 



Tsiganis et al. 



20051 ). This "numerical artifact" is useful in indicating that the timescales of 
clearing in the strongly unstable zones is < 1 My. Below we discuss in detail 
the loss of asteroids from the more stable regions of the main belt. 

3.1. Historical population of large asteroids. 

We used the test particle loss history of Sim 1 to estimate the loss history 
of large asteroids from the main belt and the large asteroid impact rate on the 
terrestrial planets. To do this we scaled /, the fraction of surviving particles 
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Figure 2: Loss history of test particles in the main asteroid belt region of the solar system 
from both Sim 1 (5760 particles for 4 Gy) and Sim 2 (115200 particles for 1.1 Gy). The 
left-hand axis is the fraction of the original test particle population that have survived the 
simulation at a given time. The right-hand axis is the estimated number of large asteroids 
in the asteroid belt, and is computed by equating the fraction remaining at t ~ A Gy with 
the number of observed H < 10.8 asteroids. The observational sample used is the 931 
asteroids with H < 10.8 excluding members of collisional families. 
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Figure 3: Removal statistics of the test particles in Sim 2 (115200 particles) as function 
of their initial orbital elements. In the upper panels, the points are colored to indicate 
their lifetime in the simulation, with red being the shortest-lived particles and blue being 
the longest-lived particles (particles surviving at the end of the simulation were removed 
for clarity). In the lower panels, the points are colored to indicate the direction in which 
they are lost: red indicates loss due to either a close encounter with Mars or removal at 
the inner boundary at 1 AU, blue indicates loss due either a close encounter with a giant 
planet or removal at the outer boundary at 100 AU, and black indicates particles that 
survived the entire 1.1 Gy simulation. 
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in Sim 1 at t = 4 Gy, to the number of large asteroids in the current main belt. 
For this purpose, we define "large asteroid" as an asteroid with D > 30 km. 
For most asteroids, size is not as well determined as absolute magnitude. 
If the asteroid visual albedo, p^, is known, the absol ute magnitude can be 
conve rted into a diameter with the following formula (IFowler and Chillemil . 



1992|)- 

p = l^lO-"/^ (1) 

Because asteroids can have a range of albedos, converting from brightness to 
diameter is fraught with uncertainty in the absence of albedo measurements. 
For simplicity, we adopt a single albedo, py = 0.09, which is approximately 



repre sentative over the size range of objects considered here (IBottke et al 



20051 ). In subsequent analysis we will use absolute magnitude as a proxy for 



size. 

Using Eq. ([1]) and our assumption of albedo, an asteroid of diameter D = 
30 km has an absolute magnitude H = 10.8. The main belt is observationally 



comp lete for asteroid absolute magnitudes as faint as = 13 (jJedicke et al. 



2OO2I ). While most H < 10.8 asteroids have existed relatively unchanged 



over the last 4 Gy, a few breakup events have created some large fragments 



over this timespan. For example, t 



lere are five members of the Vesta fam- 



ily with H < 10.8 flNesvorny et al. 



20061 ). Collisional fragments produced 



over the last 4 Gy can "contaminate" the observed H < 10.8 asteroid pop- 
ulation, and collisional breakup events have also disrupted some primordial 
H < 10.8 asteroids. These collisional processes complicate the estimate of 
the dynamical loss history of large asteroids over the age of the solar sys- 
tem. Happily, most collisional family members with H < 10.8 have been 
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identified (^Nesvorny et al. 



20061 ). and can be removed to fu rther refine the 



observation sample. The database of flNesvorny et al. 



20061 ) was also made 



with some attempt to remove interlopers that have similar dynamical prop- 
erties as a family, but a different spectral classification that indica tes they 



20051) 



are not members of the coUisional family flMothe-Diniz et al. 

The observational data set we used was the 1137 asteroids with H < 
10.8 obtained from the AstDys online data se rvice ( 



20031 ). Using the family classification system of 



Knezevic and Milanil . 



Nesvorny et al. 



fl2006h . 206 



(18%) of these asteroids are identified members of collisional families. We 
eliminated collisional family members and used the remaining sample of 931 
asteroids for our scaling. Implementing this normalization, Fig. [2] shows the 
loss history of the main belt asteroids with the population scale on the right- 
hand axis. Because the dynamical depletion of asteroids from the main belt 
is approximately logarithmic, a roughly equal amount of depletion occurred 
in the time interval 10-200 My as in 0.2-4 Gy. We find that the asteroid belt 
at t = 200 My would have had 28% more large asteroids than today, and the 
asteroid belt at t = 10 My would have had 64% more large asteroids than 
today. Our calculation indicates that ~ 2200 large asteroids [H < 10.8) may 
have been lost from the main asteroid belt by dynamical erosion since the 
current dynamical structure was established, but ~ 1600 of those asteroids 
would have been lost within the first 10 My. 
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3.2. Non-uniform pattern of depletion of asteroids 

Fig. m compares the results of Sim 1 with the observational sample [H < 
10.8 asteroids, excluding coUisional family members) ^ The proper semimajor 
axes of the surviving particles from Sim 1 at the end of the 4 Gy integratio n 
were computed using the public domain OrbitE^lcode f Knezevic et al. . 2002 ). 
The bin size of 0.015 A U was chosen using the hist o gram bin size optimization 
method described by Ishimazaki and Shinomot"o ( 2007 ) (See Appendix A). 
Fig. IDd is the ratio of the data sets. We find that the observed asteroid belt 
is overall more depleted than the dynamical erosion of an initially uniform 
population can account for, and there is a particular pattern in the excess de- 
pletion: there is enhanced depletion just exterior to the major Kirkwood gaps 
associated with the 5:2, 7:3, and 2:1 mean motion resonances (MMRs) with 
Jupiter (the regions spanning 2.81-3.11 AU and 3.34-3.47 AU in Fig. H^); the 
regions just interior to the 5:2 and the 2:1 resonances do not show significant 
depletion (the regions spanning 2.72-2.81 AU and 3.11-3.23 AU in Fig. H^), 
bu t the inner belt region ( spann ing 2.21-2.72 AU) shows excess depletion. 



Minton and Malhotral (120091 ) showed that the observed pattern of excess 



depletion is consistent with the effects of the sweeping of resonances during 
the migration of the outer giant planets, most importantly the migration of 



^Fig. 5^ is similar to Fig. la of iMinton and Malhotral (|2009l ). but with our sample of 
931 asteroids with H < 10.8 that arc not members of collisional families. The results we 



report in this section are similar to those of lMinton and Malhotral (|2009l ) ; because they are 
based on simulations with much larger number of particles, their statistical significance is 
improved. 



Found at: http://hamilton.dm.unipi.it/astdys/ 
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Jupiter and Saturn. There is evidence in the outer solar system that the giant 
planets - Jupiter, Saturn, Uranus and Neptune - did not form where we find 
them today. The orbit of Pluto and other Kuiper Belt Objects (KBOs) that 
are trapped in mean motion resonances with Neptune can be explained by the 
outward migration of Neptune due to interactions with a more m assive pri- 



mordial pla netesimal disk in the outer regions of the solar system (IMalhotral . 



1993 



19951 ). The exchange of angular momentum between planetesimals and 



the four giant planets caused the orbital migration of the giant planets until 



the outer planetesimal disk was 



planets in their present orb its (IFernandez and Ip 



1999 



depleted of most o 



its mass, leaving the giant 



1984 



Hahn and Malhotra, 



Tsiganis et al 



20051 ). As Jupiter and Saturn migrated, the locations 



of mean motion and secular resonances swept across the asteroid belt, excit- 
ing asteroids into terrestrial planet-crossing orbits, thereby greatly depleting 
the asteroid belt population and perhaps also causing a late heavy bom- 



bardment in the inner so 



2001 



Gomes et al. 



2005 



ar system ( 



Liou and Malhotra 



Strom et al. 



20051). 



1997 



Levison et al. 



We identified six zones of excess depletion; these are labeled I-IV in 
Fig.|4)D. Zone I would have experienced depletion primarily due to the sweep- 
ing uq resonance, with some contribution possibly from the 3:1 MMR. Zones 
II and IV are the zones that lie on the sunward sides of the 5:2 and 2:1 
resonances, respectively, and are hypothesized to have experienced the least 
amount of depletion due to sweeping mean m otion resonances (MMRs) an d 



Minton and Malhotral fl2009f ). 



secular resonances under the interpretation of 
Zones III, V, and VI are on the Jupiter-facing sides of the 5:2, 7:3, and 2:1 
MMRs, respectively, and are hypothesized to have experienced depletion due 
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Figure 4: Comparison between the semimajor axis distribution of Sim 1 test particles 
and the sample of observed main belt asteroids, a) The observed main belt asteroid 
distribution, Nobs, for H < 10.8 asteroids (shaded) and the surviving particles of Sim 1, 
Nsim, (solid), b) Ratio of the data sets. 



to the sweeping of these resonances. The average ratio between the model 
and observed population per 0.015 AU bin in each zone is quantified in Fig. [51 

3. 3. Empirical models of population decay 

The loss rate of small bodies fro m various regions of the so lar system has 



been studied by several authors (see 



Dobrovolskis et a. 



hensiye review of the recent literature on the subject). 



20071. for a compre- 



Holman and Wisdom 



(I1993I ) found that the decay of a population of numerically integrated test 
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Zone 



Figure 5: The average depletion, (Nobs / Nsim) , in each of the six zones identified in Fig. [lb; 
the average is taken over the 0.015 AU bins present in each zone. 
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particles on initially circular, coplanar orbits distributed throughout the 
outer solar system was asymptotically logarithmic, that is, h oc where n 
is the number of test partic les remaining in the simulation at a given time 

(120071 ) showed that, for many small body populations. 



t 



Dobrovolskis et al. 



loss is described as a stretched exponential decay, given by the Kohlrausch 
formula, 

/ = exp (- [t/t,f) , (2) 
where / is the f raction remaining o f the initial population (/ = n{t)/Ntot)- 



In the cases that 



Dobrovolskis et al 



studied, namely loss rates of small body 
populations orbiting giant planets, they found that (3 ~ 0.3. For reference we 
note that a value of /3 = 1/2 is expected for a diffusion-dominated process for 
the removal of particles; in this plot of log of the number of remaining 
particles vs. square root time would be a straight line. 



In Fig. O we adopt a similar plot style as in 



Dobrovolskis et al. 



mm 



(their Fig. 1) as a way of evaluating various empirical decay law s for the 



Dobrovolskis et al 



stretched 



results of Sim 1. Unlike the cases explored by 
exponential decay with (3 ~ 0.3-0.5 is a very poor model for the asteroid belt. 
Fig. [6] suggests either logarithmic or power law decay would be better models 
of particle decay from this simulation for t > 10^ y. Using a logarithmic decay 
law of the form: 

f = A-Blnit/lY), (3) 
and a power law decay of the form: 



/ = C{t/1 y)-^. 



(4) 



where A, B, C, and D are positive and dimensionless constants, we can look 
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Figure 6: Comparison of empirical decay laws for the main asteroid belt region from Sim 1 
(5760 particles for 4 Gy). The result from Sim 1 is plotted in five different ways. Bottom 
curve (black): fraction / of particles surviving (left-hand scale) vs. time (bottom scale). 
Next-to-lowest curve (green): log/ (right-hand scale) vs. elapsed time (bottom scale). 
Middle curve (red): log / (right-hand scale) vs. \Jt (top-scale). Next-to-uppermost curve 
(blue): / (left-hand scale) vs. logt (interior scale). Top curve (yellow): log/ (right-hand 
scale) vs. logt (interior scale). Only the yellow and blue curves resemble straight lines in 
this format, and only for t > 10^ y. 

for the logarithmic, power law, and stretched exponential functions that best 
fit the decay history of Sim 1 for t > 1 My. The best fit parameters for each 
of these decay laws are listed in Table [TJ Note that the best fit exponent D 
for the power law decay is close enough to zero that it is not very different 
than the logarithmic decay over the range of timescales considered here. 

The difference between various empirical decay laws and the simula- 
tion output from Sim 1 is shown in Fig. [71 The format of Fig. [7] is sim- 



ilar to that of Fig. 2 of 



Dobrovolskis et al. 



(120071 ). but here the y-axis is 
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Table 1: Best fit decay laws for Sim 1 (5760 particles for 4 Gy) 



Decay law Parameters Valid range (yr) 



Stretched exponential (Eq. [2]) 


log 


to = 8.6986 ± 0.0057 


t > 10^ 




/? = 


= 0.1075 ±0.0004 




Logarithmic (Eq. [3]) 


A -- 


= 1.1230 ± 0.0020 


t > 10^ 




B -- 


= 0.0377 ±0.0001 




Power law (Eq. H]) 


C -- 


= 1.9556 ±0.0027 


t > 10^ 




D -- 


= 0.0834 ±0.0001 




Piecewise logarithmic (Eq. [5]) 


Ai 


= 1.3333 ±0.0006 


10^ ° < t < 10^-2 




Bi 


= 0.05130 ±0.00004 






A2 


= ± {B2 - 5i) ■ 7.2 


10^-2 < t < 108-3 




B2 


= 0.02695 ±0.00011 






A-i 


= A2 + {B3 - B2) ■ 8.3 


10*^-3 < t < lo^-i 




Bs 


= 0.02695 ±0.00011 






A, 


= A3 ±(54-^3) -9.1 


10^-^ < t < 10^-^ 




B4 


= 0.03079 ±0.00018 





A log I In /I = log|ln/sjm| — log | ln//jt|, where the subscripts sim and fit 
refer to the simulation data and best fit model, respectively. A perfect fit 
would plot as a straight line with A log | In /| = 0. The power law function is 
a good fit, but with a small exponent D, which makes it practically similar to 
a logarithmic decay. The best fit stretched exponential value of (3 obtained 
here is much s maller than the value of ^ 0.3-0.5 obtained by many of the 



cases shown by 



Dobrovolskis et al 



(I2OO7I ). This may indicate that classical 



20 



diffusion does not dominate the loss of asteroids from tlie main belt. 

From Fig. [7] there is no clear preference for one or the other functions for 
the decay model. However, with some experimentation, we found that an 
improved fit can be obtained by considering a piecewise logarithmic decay of 
the form: 

f^ = A- Bi \n{t/l y), U<t< U+i, (5) 

where Ai and Bi are positive coefficients. A physical justification for a piece- 
wise logarithmic decay law is outlined in the following argument. If the 
region under study were divided into smaller subregions, and the loss of par- 
ticles from each of those subregions follows a logarithmic decay law, then 
the linear combination of the decay laws for all subregions is the decay law 
for the total ensemble of particles, and is itself logarithmic. However, if any 
subregion completely empties of particles, then that region remains empty 
(it cannot have negative particles), and so the decay law of that region no 
longer contributes to the decay law for the total ensemble of particles; the 
decay of the ensemble then undergoes an abrupt change in slope. A piecewise 
logarithmic decay law for an ensemble of particles originating from the main 
asteroid belt region implies that the intrinsic loss rate from the asteroid belt 
is best described as n oc but with different proportionality constants for 
different regions inside the belt. 

We found that the model that minimized A log | In /| for t > 10^ y is a 
four component piecewise logarithmic decay law with slope changes at lO'''^ y, 
10^ '^ y, and 10^'^ y. We performed a least-squares fit to the loss history of 
Sim 1, fitting it to the four component piecewise function given by Eq. ([5]); 
the best fit parameters are given in Table [H The residuals for the piecewise 
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logarithmic decay law are much reduced, compared to the other empirical 
models considered, as shown in Fig. [71 
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Figure 7: Differences between tlie best fit decay models and the loss history of Sim 1. 
Here the y-axis is A log | ln/| = log | Infsiml — log | \nffit\, where the subscripts sim and 
fit refer to the simulation data and best fit loss function, respectively. 



4. Large asteroid impacts on the terrestrial planets 

Although the impact history of the terrestrial planets is numerically dom- 
inated by small impactors, D < 10 km, the larger but infrequent impactors 
are also of great interest as they cause the more dramatic geological and 
environmental consequences. The dynamical origins of the latter have been 
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less well studied because of the unavoidable small number statistics issues 
with them. Unlike the Yarkovsky effect, which is primarily responsible for 
populating the NEA population with D < 10 km asteroids, the dynami- 
cal chaos in the asteroid belt is a size independent process and dynamical 
erosion is the primary loss mechanism for large asteroids. iMigliorini et al 



(119981 ) investigated how D > 5 km asteroids from the main belt become ter- 
restrial planet-crossing orbits, and found that weak resonance in the inner 
solar system were likely responsible for populating the NEOs. We used our 
simulations to quantify the impact rates for the larger asteroidal impactors. 

As shown in §3.31 the long term dynamical loss of asteroids from the 
main belt is nearly logarithmic in time. The fate of any particular aster- 
oid (the probability that it will impact a particular planet, the Sun, or be 



ejected from the solar system) is strongly 



in the main belt (IMorbidelli and Gladman 



dependent on its source region 



1998 



Bottke et al. 



2000). For 



example, particles originating from the region near the secular resonance 
at the inn er edge of the main belt have a ^ 1-3% chance of impacting 



the Earth f Morbidelli and Gladman 



1998 



Ito and Malhotra 



20061 ). whereas 



objects originating f urther out in the aster oid belt have much lower Earth 



impact probabilities (iGladman et al 



19971 ). Our simulations indicate that 



large asteroidal impactors that enter the inner solar system may originate 
throughout the main belt, so we needed to compute the overall impact prob- 
abilities for impactors originating by dynamical chaos from the main belt as 
a whole. We do this by means of an additional numerical simulation that 
yields the terrestrial planet impact statistics for those particles of Sim 1 that 
were 'lost' to the inner solar system. We then combine the impact proba- 
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bilities with the 4 Gy loss history of large asteroids (Fig. [2]) to estimate the 
number of large impacts onto the terrestrial planets. The details of these 
calculations are described below. 

4.I. Impact probabilities 

In our initial simulations we did not follow any of the particles all the 
way to impact with any of the planets; particles that entered the inner solar 
system were stopped either at the Hill sphere of Mars or at an inner bound- 
ary of 1 AU heliocentric distance. To compute the impact probabilities for 
the terrestrial planets, we performed an additional simulation using the re- 
sults of Sim 1. Only inward-going particles from Sim 1 were considered, as 
outward-going ones are overwhelmingly likely to collide with or be ejected 
from the solar system by Jupiter, and as Figs. [3] and illustrate, there 
are few observed asteroids beyond 3.4 AU, where outward-going asteroids 
dominate. We first identified a set of "late" particles from Sim 1 that were 
removed after 1 Gy at an inner boundary (either at the cutoff at 1 AU or by 
crossing the Hill sphere of Mars). The Kirkwood gaps are mostly emptied 
of asteroids by 1 Gy, as shown in Fig. [3l Fig. [8] shows the initial semimajor 
axes, eccentricities, and inclinations of these 136 particles. Note from this 
figure that asteroids lost due to dynamical erosion after 1 Gy can come from 
nearly anywhere in the main belt. 

First we captured the positions and velocities of the 136 late particles 
from Sim 1 at the time of their removal. We cloned each of the late parti- 
cles 128 times, such that Tdone = (1 + Sr)roriginal and Vcion.e(l + S^)Voriginal, 

where \6r\, < 0.001. The resulting 17408 particles were integrated using 
the MERCURY integrator with its hybrid symplectic algorithm capable of 
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Figure 8: Distribution of initial orbital elements of the "late" particles from Sim 1. These 
particles left the asteroid belt after 1 Gy. Only particles that were removed at an inward- 
going boundary are shown here, that is they were removed from Sim 1 either by crossing 
the inner barrier at 1 AU or crossing the Hill sphere of Mars. 
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following particles through close encounters with planets (IChambersl . 119991 ). 
In this simulation, which we designate Sim 3, all eight major planets were 
included, the nominal integration step size was 2 days, and an accuracy pa- 
rameter of 10~^^ was chosen. The particles were removed if they approached 
within the physical radius of the Sun or a planet, or if they passed beyond 
100 AU. The simulation was run for 200 My. 

For every particle that was removed in Sim 3 (either by impact or escape) 
we determined from which of the 136 source particles from Sim 1 it was 
cloned. We used the initial semimajor axis of a given source particle and 
placed it into a bin of 0.015 AU in width, which we call the source bin. 
We then weighted the removal event by a factor equal to the ratio of the 
abundance of observed H < 10.8 asteroids to the abundance of particles 
at the end of Sim 1 in the source bin. The weighting factor, which also 
quantifies the relative amounts of depletion throughout the asteroid belt, is 
shown in Fig. IDd. This weighting accounts for the differences in the orbital 
distribution of our simulated asteroid belt and the observed large asteroids. 
The raw impact statistics as well as the probabilities weighted based on 
the distribution of observed asteroids are tallied in Table [2j These give the 
probability that an asteroid originating in the main belt, and which becomes 
a terrestrial planet-crosser by dynamical chaos, will impact a planet, the Sun, 
or be ejected from the solar system. The unweighted and weighted terrestrial 
impact probabilities for the terrestrial planets differ by < 10%, and both give 
an impact probability onto the Earth of ~ 0.3%. 
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Table 2: Impact probabilities for the terrestrial planets for Sim 3. Ejected (OSS) refers to 
particles that either crossed the outer barrier at 100 AU or encountered the Hill sphere of 
a giant planet. 



Fate 


Number 


% 


Weighted % 


Ejected (OSS) 


13862 


78.6 


70.0 


Survived 


269 


1.55 


13.7 


Sun 


3111 


17.9 


15.3 


Mercury 


10 


0.057 


0.061 


Venus 


74 


0.425 


0.396 


Earth 


56 


0.322 


0.306 


Mars 


26 


0.149 


0.140 



4-2. Flux of large (D>30 km) impactors on the terrestrial planets 

We used the weighted impact probabihties shown in Table [2] and our 
model of the loss history of asteroids from Sim 1 to estimate the number of 

> 30 km impacts onto the terrestrial planets since the end of the LHB. 
Here we make the assumption that t = in our model is about 4 Gy ago, 
roughly the post-LHB era. Because the loss rate of asteroids is approximately 
logarithmic, the loss rate is much higher at early times than later ones. The 
early time is likely to have coincided with the tail end of the LHB itself; 
separating out the component of the impact flux that is due to the LHB 
rather than to dynamical erosion is problematic. We therefore consider only 
the dynamical loss at t > 100 My in our model as part of the post-LHB epoch. 
With these assumptions, our model finds that since the end of the LHB, the 
Earth has experienced ~ 1 impact of a D > 30 km asteroid. Venus, with 
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its slightly higher impact probability than Earth, should have experienced 
~ 1.3 impacts of this size since the end of the LHB. Our model also suggests 
that Mars and Mercury have had ~ 0.6 and ~ 0.1 impacts of > 30 km 
asteroids, respectively. These small numbers are consistent with there having 
been no impacts of D > 30 km asteroids on the terrestrial planets since the 
end of the LHB. 



4-3. Comparison with record of large impact craters on the terrestrial planets 
Known large impact basins on the terrestrial planets are generally of 



ages confined t o the first 800 



Bombardment (IHartmann 



1965 



My of s olar s ystem history, the Late Heavy 



Ryder 



2OO2I ). Excluding those, we consider 



only the post-LHB large craters on Earth and Venus. 

The three largest known impact structures on Earth are Vredefort, Sud- 
bury, and Chicxulub crater s, each with f i nal cr a ter diameters Dr < 300 krr i 



and ages less than ~ 2 Gy (IGrieve et al. 



20081). 



Turtle and Pierazzd (119981 ) 



argue that Vredefort crater has a diameter of Dc < 200 km, which would 
make the final diameters of all of the largest known impact craters on Earth 
Dc ~ 130-200 km. At least one of these large impact events has been associ- 
ated with a mass extinction. The Chixculub crater, estimated to have been 
created by the impact of a ~ 10 k m obiect, is associated with the term! 



nal Cretaceous mass extinction event (Alvarez et al. 



19801: 



Hildebrand et al. 



I991I ). Estimates from impact risk hazard assessments in the literature sug- 



gest that Chixcul ub- sized impact events happen on Earth on the order of 



1994 ) 



once every 10* y (1 Chapman and Morrisonl . 

The impact cratering record of Venus is unique in the solar system. Over 
98% of the surface of Venus was mapped using synthetic aperture radar by the 
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Magellan spacecraft (iTanaka et al.l . 119971 ) . The observed craters on Venus ap- 
pear mostly pristine and are randomly distributed across the planet's surface, 
which has been taken as evidence fo r a short-lived, global resurfacing eyent o n 



the planet within the last ~ 1 Gy fjPhillips et al 



1992 



Strom et al 



19941 ) 



Venus has four craters with > 150 km; Mead crater with = 270 km, 

Isabella crater with Dc = 1 75 km, Meitner with D r = 149 km, and Klen- 



ova with Dc = 141 kmo 



Shoemaker et al. 



(jl99ll ) estimated the surface 



age of Venus to be ~ 200-500 My using the total abundance of Venus 
craters and a flux of impactors based on the known abundance of Venus- 
crossing asteroids and on mode ls of their impact probabilities. More recently, 
Korycansky and Zahnld (120051 ) used similar techniques (as well as an atmo- 
spheric screening model for small i mpactors) and estimated the surface age 
of Venus to be 730 ± 220 My old. Iphillips et al.l Jl992l ) used four methods 
to determine the age of Venus' surface, three which used models of observed 
Venus-crossing asteroids and one that used the observed abundance of craters 
on the lunar mare as a calibration. All four methods resulted in surface ages 
between 400-800 My. 

Each of the techniques described above for estimating surface ages based 
on the abundance of observed craters has its shortcomings. Calculating a 
surface age using the observed population of NEAs makes the assumption 
that the current population of NEAs is typical for the entire post-LHB history 
of the inner solar system, including both the number of near earth asteroids 
and their computed impact probabilities. Calculating the surface age based 



^From the USGS/University of Arizona Database of Venus Impact Craters at 



http: / / astrogeology.usgs.gov/Projects/VenusInipactCraters/ 
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on the abundance of craters on the lunar mare makes the assumption that 
the crater production rate in the inner solar system has been approximately 
constant over t he last ^ 3.2 Gy, which is the age by when most lunar mare 
were produced JbVSP. 198i[ ). 

Large impactors produce craters with a complex morphology. To deter- 
mine the sizes of impactors which created the largest post-LHB terrestrial 
impact craters, we used Pi scaling relationships to estimate the transient 
crater diameter as a function of projectile and target parameters; we then 
applied a second scaling relationship between transient crater diameter and 
final crater d i amete r. The particular form of Pi scaling used here is given by 
(120051 ) for impacts into competent rock: 



Collins et al 



D 



tc 



1.161 



1/3 



^0.78^0.44^-.22gj^ 



(6) 



where pi and pt are the densities of the impactor and target in kg m""^, D 
is the impactor diameter in m, Vi is the impactor velocity in m s"^, g is 
the acceleration of gravity in m s~^ and 9 is the impact angle. We used the 
relationship b etween final crater diarneter, and transient crater diameter, 
Dtci given by lMcKinnon and SchenkI (119851 ): 

^1.13 



= 1.17 



'tc 



D 



0.13 ' 



(7) 



where is the diameter at which the transition between simple and com- 
plex crater morphology occurs. The transition diameter, D^ , is inversely 
proportional to the surface gravity of the target (jMeloshi . Il989l ). and can be 
computed based on the nominal value for the Moon of D^^ ^oon = 18 km. 

Applying Eqs. ([6]) and ([7]) to the problem of terrestrial planet cratering by 
asteroids requires some assumptions about both the impacting asteroids and 
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the targets. For simplicity we assumed that target surfaces have a density 
Pt = 3 g cm~^ and that impacts occur at the most probable impact angle of 



45° fiGilbert 



18931 ). The characteristic impact velocity is often chosen to be 



the rms im pact velocity o btained from a Monte Carlo simulation of planetary 
projectiles (IBVSP. 198ll ). However, using the rms impact velocity may lead 
to misleadingly high estimates of the i mpact velocity because the impact ve 



locity distributions are not Gaussian teottke et al. 



1994 



Ito and Malhotra 



20091 ). Therefore the media n veloci t y may be m ore appropriate estimate of 



a "typical" impact velocity (IChyba 



1990 



19911 ) 



We made our own estimates of the impact velocities onto the planets 
using the results of Sim 3. The small number of impacts computed in Sim 3 
makes determining impact velocity distributions difficult. We improved the 
statistics by using the far more numerous close encounters between planets 
and test particles to calculate mutual encounter velocities. For every close 
encounter in Sim 3 we recorded the closest-approach distance and mutual 
velocity between particles and planets. We only considered particles that had 
a closest-approach distance within the Hill sphere of a planet. We estimated 
impact velocities using the vis viva integral given by: 



1 , GM„ 1 . GM„ 
—V = —V 

9 r n enc ' 



(8) 



where Vimp and Tp are the estimated impact velocity and the radius of the 
planet, and Venc and Tenc are the mutual encounter velocity and closest ap- 
proach distance. Some particles encountered a planet multiple times, which 
skews the impact velocity estimates, so we only considered unique encounters 
between a particular particle and a planet. If a particle encountered a planet 
multiple times it was treated as a single event and we only used the impact 
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velocity of the first encounter. The median, mean, and RMS impact veloc- 
ities for each of the planets is shown in Table O The velocity distributions 
are shown in Fig. [91 



i: Estimated impact velocities for cfose encounter events in 


Planet 


Encounters 


Velocity (km 








Median 


Mean 


RMS 


Mercury 


3885 


38.1 


40.5 


43.3 


Venus 


6974 


23.4 


25.9 


27.5 


Earth 


3522 


18.9 


20.3 


21.1 


Mars 


1205 


12.4 


13.1 


13.8 



Using the median impact velocities and our assumptions about asteroid 
density (pj = 1.5-3 g cm~^), we calculated the sizes of the impactors that pro- 
duced the largest post-LHB craters in the inner solar system, using Eqs. 
and ([7]). The "big three" terrestrial impact craters are consistent with aster- 
oidal impactors of diameter D ~ 7-14 km. The estimated projectile size for 
Mead crater, the largest impact crater on Venus, is D ~ 12-16 km. Varying 
the assumptions used in the Pi scaling. Mead and Isabella craters are fully 
consistent with an impact of a D > 10 km asteroid, but Meitner and Klenova 
craters could be consistent with smaller asteroidal impacts. 

In summary, there are no known impact structures of ages < 3.8 Gy (post- 
LHB) attributed to projectiles with D >30 km. This is consistent with the 
theoretical estimate above. The small number of observed large impacts are 
consistent with the impact of D > 10 km objects. Earth has been impacted 
by at least 3 objects that are consistent with D > 10 km asteroids. During 
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Figure 9: Impact velocity distributions of asteroids on the terrestrial planets. 
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the Phanerozoic eon (the last 545 My), only one impact crater, Chixculub, 
has been discovered on Earth that is consistent with a, D > 10 km asteroid. 
Venus has been impacted by 2-4 D > 10 km objects in the last Gy, based 
on the number of observed large craters and its estimated surface age. 

4.4- Flux of D>10 km impactors on the terrestrial planets 

In our discussion of the dynamical erosion of the main asteroid belt, we 
have confined ourselves to D > 30 km primordial asteroids because non- 
gravitational and collisional effects are negligible for this population. How- 
ever, considering that the largest craters on the terrestrial planets correspond 
to impactors D ~ 10 km, somewhat smaller than 30 km, we are motivated 
to consider the > 10 km population of the main belt. In the size range 
D = 10-30 km the effects of collisions and non-gravitational forces are not 
negligible, but they do not dominate that of dynamical chaos on the semima- 
jor axis mobility of main belt asteroids. Therefore we extend our dynamical 
calculations to > 10 km asteroids and we compare the results with the 
terrestrial planet impact crater record, with the caveat that our results are 
only a rough estimate of the actual impact flux. 

In order to turn the asteroid loss rate into an impact flux for a given 
crater size, the results of Sim 1, previously normalized to the abundance of 
H < 10.8 asteroids (Fig. [2]), were scaled to the fainter asteroids (if < 13.2). 
Applying = 0.09 for the geometric albedo, H < 13.2 corresponds to 
-D > 10 km. The size distribution of the large asteroids of the main belt 



(Bottke et al.. 


2005; 


Strom et al.. 


2005) 



took the H < 10.8 {D > 30 km) loss rate shown in Fig. [2] and scaled it to 
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H < 13.2 [D > 10 km) assuming the asteroid cumulative size distribution in 
this size range is a simple power law of the form: 

NyD = kD-\ (9) 

The debiased main belt asteroid size frequency distribution determined by 



Bottke et al. 



(120051 ) gives h = 2.3 in the size range 10 km < D < 30 km. 
We used our estimate of the loss of -D > 10 km asteroids from the main 
belt due to dynamical diffusion, along with our estimate of the impact prob- 
abilities shown in Table [21 to determine the impact flux of > 10 km 
asteroids, A^>io kmj on the terrestrial planets. The impact flux is based on 
the four-component piecewise decay law given by Eq. ([5]) with parameters 
given in Table [H which gives the fraction remaining as a function of time. 
First the derivative of Eq.([S]) is taken, yielding f{t). In order to convert 
from a fraction rate, /(t), to the flux of > 10 km impacts, A^>io km(^); we 
multiplied f{t) by the coefficient, C>io km, defined as: 

931 /lO^ "^-^ 



The first component of the coefficient is the constant, 931//(4 Gy), which 
normalizes the fraction remaining such that the total number of asteroids 
remaining at t = 4 Gy is 931; the latter is the total number of H < 10.8 
(-D > 30 km) asteroids in the observational sample. The next component 
is (10/30)"^'^, which scales the results to > 10 km, as given by Eq. Q. 
Finally, the coefficient was multiplied by the weighted probability Pi that 
terrestrial planet-crossing asteroids impact a given planet. For the Earth, 
the impact probability is PEarth = 0.003 as given in Tabled By multiplying 
f{t) by the coefficient C>io km we obtained the flux of impacts by -D > 
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10 km asteroids on the Earth as a function of time, A^>io km- By further 



mu 

by 



■tiplying by 1/23, 



using the ratio of 23:1 Earth:Moon impacts calculated 



Ito and Malhotral fl2006l ). we also estimated the lunar flux. The results 
are shown in Fig. [TOk . showing our estimate of the impact flux of £) > 10 km 
asteroids on the Earth (left-hand axis) and Moon (right-hand axis) since 
1 My after the establishment of the current dynamical architecture of the 
main asteroid belt. 



For co mparison, we a . 



tained by 



Neukum et al 



so p 



ot in Fig. [TUk the impact flux estimates ob- 



200ll ) from crater counting statistics. The shaded 



region of Fig. [TOk represent upper and lower bounds on the estimated post 



LHB impact flux using calibrated lunar cratering statistics ( iNeukum et al. 



200ll ). The lower bound is calculated using the number of > 200 km 



lunar impact c r aters from the Neukum Production Function (NPF, Fig. 2 of 



Neukum et al. 



( 20011 1): the rate is 7 x 10"^ craters km~^ Gy~^. The upper 
bound is calculated using the number of Dc > 140 km lunar impact craters 
from the NPF; the rate is 2 x 10"^ craters km~^ Gy~^. These rates were 
multiplied by the surface area of the Moon to obtain the number of impacts 
on the lunar surface, and then multiplied by 23 to obtain the number of 
impacts on the Earth's surface. 

We have also calculated the cumulative number of impacts on a surface 
with a given age; the result is shown in Fig. [TOb . The cumulative num- 
ber of impacts as a function of time is calculated simply as Ncumuiative = 
C>io km [fit) — /(4 Gy)]. Assuming that t = corresponds with an age of 
4 Gy ago, the surface age is simply SA = 4 Gy — t. The upper and lower 
bounds on the cumulative number of impacts calculated from the NPF are 
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shown as the shaded region, similar to Fig. [TUh . Our calculations of the im- 
pact flux show that at early times the flux was larger than the estimate given 
by the NPF, but that after t = 500 My the flux of impacts that we calculate 
is lower than that given by the NPF. The present day flux of i5 > 10 km 
impactors is currently about an order of magnitude lower than that given by 
the NPF. 

We note that the dynamical erosion rate decays with time, and the impact 
rate that we derive also decays as a function of time. Our estimated impact 
flux for large craters declines by a factor of 3 over the last 3 Gy, as shown 
in Fig. [TOl This is in some contrast with previous studies: in commonly 
used cratering chronologies, the impact flux is usually assumed to have been 
relatively constant over time since the end of the LHB. The results obtained 
for the Earth-Moon system shown in Fig. [10] can also be applied to the 
remaining terrestrial planets, using the weighted impact probabilities, pi, 
given in Table [H Our estimated decay of the large impa c t rate for Mars over 



Ouantin et al. 



(120071 ) that suggest 



the last 4 Gy agrees with the results of 
that the impact cratering rate of -De > 1 km craters on Mars has declined by 
a factor of 3 over the last 3 Gy, based on counts of craters on 56 landslides 
along the walls of Valles Marineris. Reliable estimates are lacking for the 
absolute ages of lunar surfaces with ages < 3 Gy, but what estimates do 
exist (i.e. estimated ages of Copernicus, Tycho, North Ray, and Cone craters 
from the Apollo missions) seem to be mostly consistent with a constant flux 
of impacto rs after 3 Gy ago, but w ith a possible increase in the flux after 



1 Gy ago (IStoffler and Ryderl . 



200ll ). Also, estimates of the current impact 



flux from studies of the lunar record seem to be consistent with estimates 
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made by obse rving; the modern NEO populatio n and estima t ing th eir impact 



probabihties (IMorrison et ah 



2OO2I). However, 



Culler et al 



f l2000h suggested 



that the lunar impact flux declined by a factor of 2 or 3 from ~ 3.5 Gy ago 
until it reached a low at ~ 600 My ago and then increased again, based on 
dating of lunar impact glasses in soils. This hypothesis is consistent with our 
result, that the overall impact flux onto the terrestrial planets has decreased 
by a factor of 3 since ~ 4 Gy ago due to the dynamical erosion of the asteroid 
belt. The apparent increase in the flux at ~ 600 My ago until the present may 
be the result of a few large asteroid breakup events in the inner asteroid belt, 
such as the Flora and Baptistina family-forming events, and therefore the 
modern NEO population may not be representative of the NEO population 
over the last ~ 4 Gy. 

Assuming that t = corresponds to an age of 4 Gy ago, we estimate 
that in the last 3 Gy there have been ~ 4 impacts of > 10 km asteroids 
on the Earth due to dynamical erosion, and in the last 1 Gy there has been 
only 1 impact. This is about an order of magnitude lower than estimates of 
the cratering rate of the terrestrial planets using the lunar cratering record, 
which give ~ 6-18 -D > 10 km impactors onto the Earth every 1 Gy us- 
ing the NPF and assuming a -D > 10 km object pr o duces a 150-200 km 



crater ( Neukum et al. 



2001 



Chapman and Morrisonl . Il994l ). The discrep- 



ancy between our estimates of the production of large craters and those 
based on the NPF may indicate that either (1) the current production rate 
of large craters {Dc > 150 km) is substantially overestimated by the NPF, or 
(2) the current production rate of large craters is not dominated by chaotic 
transport of large main belt asteroids. 
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Figure 10: Impact rates of Z) > 10 km asteroids on the Earth and Moon, a) Sohd hne 
is the impact flux of D > 10 km impactors per Gy onto the Earth and the Moon, from 
our calculations, b) Solid line is the cumulative number D > 10 km impactors for a given 
surface ag e, from our calculatio ns. The shaded regions indicate the Neukum Production 



Function (Neukum et al 



200lh upper and lower bounds on the production of Dc > 140- 



200 km impact craters. 
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One possibility is that cratering at this size is dominated by comets, 
rather than asteroids. The fraction of terrestrial planet impacts that are due 
to comets vs. asteroids has long been controversial, but is generally thought 
that Jupiter-family comets c ontribute few e r than 10% and Port cloud comets 



contribute fewer than 1%. flBottke et al. 



2002a 



Stokes et al. 



20031 ). Even 



periodic comet sh owers may not substant ially increase the impact contribu- 



tion from comets (iKaib and Quinnl . [200^. Therefore it seems unlikely that 



comet impacts can account for an order of magnitude more large impacts on 
the terrestrial planets than asteroids. 

Another possibility is that large asteroid breakup events (followed by 
fragment transport via the Yarkovsky effect to resonances) dominate the 
production of large craters. Our model neglects breakup events which have 
produced numerous ~ 10 km fragments over the last 4 Gy, and breakup 
events near resonances with Earth-impact probabilities higher than that of 
the intrinsic main belt may contribute to the large basi n impa ct rate. This 



Bottke et al 



(120071 ) for the ori- 



is similar to the hypothesis proposed by 
gin of the Chicxulub impactor from the Baptistina family-forming event. 
The Baptistina breakup is hypothesized to have involved two large asteroids 
{Di ~ 170 km and D2 ~ 60 km) that collided at a semimajor axis distance 
< 0.01 AU from two overlapping weak resonances that, combined, increased 
fragment ec centricities to plane t-crossing orbits with a 1.7% Earth-impact 
probability (IBottke et al.l . 120071 ). It is unclear whether such fortuitous com- 
binations of conditions occur often enough to dominate the production of 
large craters over the solar system's history since the end of the LHB. Large 
family forming breakup events do occur in the asteroid belt, and they likely 
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can increase the flux of impacts onto th e terrestrial planets (INesvorny et al. 



2002 



2006 



Koro chant seva et al 



20071 ). The Flora family breakup event in 



part icular was large and oc curred in the inner asteroid belt roughly 0.5-1 Gy 



ago flNesvorny et al. 



2OO2I ). Also, while the Yarkovsky effect is very weak for 



D > 10 km asteroids, it is not entirely negligible over the age of the solar sys- 
tem. For instance, the loss rate of D = 10 km from the inner asteroid belt is 
somewhat higher when the Yarkovsky effect is taken into account compared 



to when it is not (IBottke et al. 



2002bl ). It is doubtful whether this modest 



difference can account for a factor of ten increase in the flux of -D > 10 km 
objects in the terrestrial planet region, however additional modeling may be 
needed to conflrm this. In addition, the weakness of the Yarkovsky effect on 
large asteroids can paradoxically enhance their mobility in the inner main 
belt. Combinations of nonlinear secular resonances and weak three-body res- 
onances m ay cause astero i ds to slowly diffuse through th e middle and inner 



main belt ( Carruba et al. 



2005 



Michtchenko et al. 



2OO9I ) . Only large aster- 



oids that are only weakly affected by the Yarkovsky effect can remain inside 
these resonances long enough for them to act. 

Finally, we consider the contribution from high velocity impactors with 
-D < 10 km to the large impact crater production rate. The velocity distri- 
butions of asteroid impacts on the terrestrial planets have signiflcant high 
velocity tails, as seen in Fig. [91 Because size distributions of asteroids fol- 
low a power law with a negative index, as in Eq. smaller objects are 
more numerous than larger ones. We calculated the relative contribution of 
impacts by objects of varying sizes on the production of craters of a given 
size. Using Eqs. ([6]) and ([7]), a D = 10 km projectile striking the Earth at 
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19 km s"^ produces a crater with a final diameter Df. = 187 km, assuming 
target and projectile densities are both 3 g cm~^, and the impact occurs at a 
45° angle. We solved for the projectile size needed to produce a Dc = 187 km 
crater while varying the impact velocity. The result was used to transform 
the Earth impact velocity distribution from Fig. [Hlinto a projectile size dis- 
tribution for a fixed final crater diameter. Finally we multiplied the binned 
projectile size distribution by a binned asteroid size distribution, using a cu- 
mulative distribution as in Eq. with an index b = 2.3. The result was 
then turned into a cumulative distribution and normalized to > 10 km, 
and is plotted in Fig. [TTl This plot shows that asteroids with D < 10 km 
that impact at high velocity increase the production of large craters by no 
more than a factor of two. Therefore D < 10 km asteroids impacting at high 
velocity cannot account for the order of magnitude difference in the produc- 
tion rate of large impact craters on the Earth between our model and the 
NPF. 

On very ancient terrains, the discrepancy between total number of ac- 
cumulated craters estimated from our model compared with the constant- 
flux models used in crater chronologies is less than with younger surfaces, 
as illustrated by Fig. [TUb . Our model also cannot account for the LHB 
itself. Total accumulated craters on ancient heavily cratered terrains asso- 
ciated with_Jlie_X;HB_are at lea st 10-15 times higher than on younger ter- 



rains ( Stoffler and Ryderl . 



200 ll ). and likely even more if the surfaces reached 



equilibrium cratering. The impact flux during the LHB was at least two or- 
ders of magnitude higher than the average flux over the last 3.5 Gy, an d possi 



bly three orders of magnitude more if the LHB was a short-lived event (IRyderl . 
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Figure 11: Cumulative size distribution of projectiles that contribute to the production of 
Dc = 187 km craters on Earth. This is calculated by solving for the projectile size needed to 
produce a crater on Earth with a final diameter Dc = 187 km for the velocity distribution 
of Fig.[SJ;, and convolving the result with the asteroid size frequency distribution (Eq. 
with index b = 2.3). This cumulative distribution is normalized such that the number of 
craters produced by objects of diameter Z? > 10 km is unity. 
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20021 ). Even with the more enhanced rate of impacts at early times, and as- 
suming higher Earth impact probabihties of 3% based on estimates from the 
resonance, the impact flux due to dynamical erosion is roughly one or two 
orders of magnitude lower than th at needed to produce the heavily craterec 



terrains associated with the LHB (IStoffler and Ryderl . 
200lk 



2001 



Neukum et al. 



5. Conclusion 

The main asteroid belt has unstable zones associated with strong orbital 
resonances with Jupiter and Saturn and relatively stable zones elsewhere. In 
most of the strongly unstable zones, the timescale for asteroid removal is < 1 
Myr; one exception is the 2:1 mean motion resonance with Jupiter where 
the timescale to empty this resonance approaches a gigayear. The relatively 
stable zones also lose asteroids by means of weak dynamical chaos on long 
timescales. We have found that the dynamical loss of test particles from the 
main asteroid belt as a whole is best described as a piecewise logarithmic 
decay. A piecewise logarithmic decay law implies that the intrinsic loss rate 
from the asteroid belt decays inversely proportional to time, n oc t~^, but 
with different proportionality constants for different regions inside the belt. 
When a region with a particular decay rate empties of asteroids, the decay 
law for the entire asteroid belt undergoes a change in slope. This logarithmic 
loss of asteroids due to dynamical chaos was established very soon after the 
current dynamical architecture of the asteroid belt was established, and it 
continues to the present day with very little deviation. Dynamical chaos is 
the predominant mechanism for the loss of large asteroids, D > 10-30 km. 
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We have calculated that the asteroid belt 1 My after the establishment of 
the current dynamical architecture of the solar system had roughly twice 
its present number of large asteroids. Because their loss rate is inversely 
proportional to time, we estimate that the flux of large impactors has declined 
by a factor of 3 over the last 3 Gy. We have calculated that large asteroidal 
impactors originating across the main asteroid belt have an overall Earth 
impact probability of 0.3%, and that the number of impacts of D > 10 km 
asteroids on Earth is only ~ 1. Our result on the current impact flux of 
D > 10 km asteroids due to dynamical erosion of the asteroid belt is an 
order of magnitude less than the values adopted in the recent literature on 
crater chronology and impact hazard assessment. We have evaluated several 
possible explanations for the discrepancy and find them inadequate. Our 
results can be used to improve studies of large impacts on the terrestrial 
planets. 

Appendix A. Determining the optimal histogram bin size 

When presenting data as a histogram, the problem of what bin size to 
choose arises. In the literature, the choice of bin size is often ad hoc, but 
need not be so. In the case of the asteroid belt, we wish to analyze the 
semimajor axis distribution of asteroids in order to gain insight into its past 
dynamical history. If we choose a bin size that is too small, stochastic vari- 
ation between bins can mask important underlying variations in the orbital 
element distributions. If we choose a bin size that is too large, important 
small-scale variations become lost (i.e. the variability near narrow Kirk- 
wood gaps). In this work, we use histogram bin size optimizer developed 
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by IShimazaki and Shinomotol (120071 ) for optimizing time-series data with 
variability that obeys Poisson statistics. This method is intuitive and easy to 
implement, and can be generalized to the problem of the number distribution 
of asteroids as a function of semimajor axis. 

The optimal bin size is found by minimizing a cost function defined as: 

C(A) = H^, (11) 

where the data have been divided into N bins of size A, /c is the mean of the 
number of asteroids per unit bin, and u is the variance. The mean is 

1 ^ 

^ = ]vE^- (12) 



1 ^ 

-rrT^^h-kr. (13) 



and the variance is 

N 

N 

The method works by assuming that the number of occurrences (or in our 
case, the number of asteroids), /c, in each bin obeys a Poisson distribution 
such that the variance of k is equal to the mean. Fig. [12] shows the C(A) 
for distribution in proper semimajor axis of our reduced set oi H < 10.8 
asteroids described in ^ The optimal bin size is ~ 0.015 AU for this data 
set. The results are similar for the surviving particles of Sim 1. 
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